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We report calculations for electronic ground states of parabolically confined quantum dots for up to 30 elec- 
trons based on the quantum Monte Carlo method. Effects of the electron-electron interaction and the response 
to a magnetic field are exposed. The wavefunctions and the ground state energies are compared with purely 
classical calculations performed with a comprehensive Molecular Dynamics code. For the chosen well parame- 
ters a close correspondence in the overall shape of electron density distribution is found even for small number 
of electrons, while the detailed radial distributions show the effects of Pauli principle in the quantal case. 
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I. INTRODUCTION 



The experimental and technological advances in production 
and use of two-dimensional quantum dot structures provide 
also new challenges and possibilities for theoretical studies 
reaching into surprisingly wide areas of physics 0] . For ex- 
ample, techniques for monitoring of controlled single-electron 



coherent transport through double dot systems has recently 
been realized 01 . Studies of structural properties of few- 00] 
and many-electron systems 0] are in this perspective 
of crucial importance for experimental progress. In particular 
coherent and quasi-coherent dynamics of quantum dots with 
a large number of electrons pose both theoretical and exper- 
imental challenges. 'Large' is here used in the sense 'larger 
than two', but still too few to justify application of mean field 
approaches |9]. Coherent dynamics and the interaction with 
a classical environment leading to decoherence are so far a 
highly unexplored area of the border between quantum and 
classical physics. 

Shell-filling effects in quantum dots were demonstrated as 
oscillations in the addition energy spectra with increasing 
number of electrons in the dot almost ten years ago 
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few years earlier the ground state energies as function of an 
external magnetic field were measured full . These findings 
sparked off an intense activity on the applications of well- 
known quantum mechanical structure theories from atomic 
and nuclear physics to the new area of quantum dots. The 
Hartree-Fock (HF) method, Linear Combinations of Atomic 
Orbitals (LCAO) methods, Density Functional Theory (DFT) 
and the quantum Monte-Carlo (QMC) methods are here the 
most frequently applied (cf. e.g. (3) working tools which 
have been used and further developed in quantum dot studies. 
In structure theory it is well known that LCAO combined with 
excact diagonalization rapidly becomes computationally im- 
practical with increasing number of electrons. HF, DFT and 
QMC theories are in general much more applicable with in- 
creasing particle size. The DFT method method has success- 
fully reproduced experimental addition energy spectra fl2ll . 
while similar success of the HF method has not yet been re- 



ported. The HF method as well as mean field models do how- 
ever give rise to oscillatory addition energy spectra ffl as 
well. Various QMC approaches have investigated the role of 
correlations, confinement geometr y an d shell effects for zero 
external magnetic field, see eg. |15l [ul [Ij, UM- The 
onset of a magnetic field introduces an additional parameter 
in the Hamiltonian and thus an implicit dependence of the N- 
electron energy levels as function of the field. It is well known 
that with increasing (weak) field the various field free energy 
levels depend differently on the field which leads to an intri- 
cate spectrum of (avoided) level crossings. It was for example 
recently shown that the character of ground state of a two- 
electron double dot changes from a S = type to a S = 1 
and back to S = type as the field strength is increased 01. 
The behavior of many-particle systems in magnetic fields has 
previously been investigated for up to 13 electrons 0,0] by 
QMC methods. 

In this paper we apply a new QMC procedure l20ll to 
study the ground state properties of parabolically confined 
two-dimensional quantum dots. We calculate addition en- 
ergy spectra for dots with 2 to 20 electrons and we dis- 
play the behavior of the ground state for increasing magnetic 
field. Further more, the energy and the probability density 
are compared with purely classical electron dynamics calcula- 
tions. These are based on the employment of non-equilibrium 
molecular dynamics methods where a thermostat slowly cools 
the system to a near-frozen state. 

In the following section we briefly describe the theoretical 
models and the computational methods. The results are dis- 
played and discussed in section III followed by concluding 
remarks in section IV. 



II. MODELS 

Our model Hamiltonian H of N identical electrons 
parabolically confined in two dimensions Ti = (xi,yi) is 
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given by 
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Here m* is the effective electron mass and e is the dielectric 
constant of the background material. Realistic experimental 
values for a two dimensional structure surrounded by a GaAs 
material is m* = 0.067ra e where m e is the electron mass 
and e = 13.5. The parameter a is empirically introduced to 
mimic a finite confinement size in the z— direction. A non- 
zero a will generally increase the mean distance between real 
confined dot electrons and thereby reduce the role of electron- 
electron interactions. The confining harmonic potential is de- 
fined by the frequency uoq which in real systems typically is 
hujQ ~ 2 — 6 meV. 

In addition the dot electrons may be exposed to a constant 
magnetic field, B, in the z— direction which correspond to an 
electromagnetic vector potential A = \B {— e x -\ — e y ). The 
effective external single particle potential can then be brought 
to the form, 
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with L z being the angular momentum operator of the z— axis. 
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The effective confining potential is given by w 2 = uOq + ^f-, 

where ujb = is the Larmor frequency. In the following 
we apply reduced atomic units and set uj = h = m* = 1. 



A. Variational Monte Carlo 

In the contruction of the single-particle basis and the trial 
wave function to be used in our variational Monte Carlo sim- 
ulations we follow closely the work of Harju et al, see for 
example Ref. 00. 

We thus assume that the electrons in the quantum dot move 
in a two-dimensional parabolic potential, meaning in turn that 
the one-body problem is similar to the classical harmonic os- 
cillator problem. In scaled units, our ansatz for the single- 
particle basis is given by 
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VV± H oc (x±%)H L M |m|/2) (r 2 )exp-(y), (3) 

with n being the shell index, m the angular momentum and 
r 2 = x 2 + y 2 . Since this is also a system with no well- 
defined center of mass motion, we need, in order to obtain 
a translationally invariant Hamiltonian to redefine the coordi- 
nates through 



x±iy=(x- x cm ) ± i(y - y cm ), 



(4) 



where x cm and y cm are the center of mass coordinates. In 
addition to the harmonic oscillator part for the single-particle 
motion we assume as discussed above that the electrons inter- 
act via repulsive Coulomb forces. 



Our ansatz for the many-body wave function consists thus 
of a Slater determinant multiplied by a two-body correlation 
function /(r^-), leading to the following form for the TV- 
particle trial wave function \I> n 
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where ipi refer to the single-particle wave functions defined in 
Eq. with the index i representing the single-particle state 
i with quantum numbers n, m and positions [xi — x cm ) ± 
i(yi — 2/cm)- The single-particle wave functions depend on x 
and y and the center of mass coordinates x cm and y cm . In 
our approach here we do not introduce variational parameters 
in the single-particle wave functions of Eq. 0, see also the 
discussions in Refs. 01. 

For the two-body correlation function we assume one of the 
simplest possible Jastrow factors, namely 



f(rij) =exp( 
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with a\ and a 2 being variational parameters. Note that these 
parameters are different for fermion pairs with parallel and 
anti-parallel spins, respectively. This leads in total to four 
variational parameters in our calculations. Note also that the 
correlation function does not depend on the center of mass 
motion. The work of Alexander and Coldwell B26I1 contains 
an extensive list of correlation functions for light atoms sys- 
tems. 

The norm of the total wave function is not needed since it 
is redundant in the Metropolis sampling. For more details on 
the Metropolis sampling and variational and diffusion Monte 
Carlo techniques see for example Refs. 0, |^ |28|] . We 
apply here a standard Metropolis sampling, as discussed in for 
example Refs. 00] where the main problem is to integrate 
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In our discussion below, we compare the energies obtained 
with the above correlated wave function with those arising 
from the non-interacting case, given by the sum over single- 
particle energies determined by the harmonic oscillator func- 
tion. These energies are given by the Fock-Darwin energy 
spectrum ll29ll 



hu(n x + riy + \m\ + 1) 
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All simulations in this work were performed with 10 T 
Monte Carlo cycles (and 5 x 10 5 thermalisation steps). There 
are two variational parameters in the Jastrow factor of Eq. ©, 
a\ = 1.0 and a 2 = 1.5 are found to give consistently op- 
timized results for all considered values of N. In Table [0 we 
show some of the simulation results for different values of N, 
with a\ and f3 as described above. 
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N 


Electron configuration 


energy 


standard deviation 


2 


TT 


3.15853 


0.000408482 


2 


TT 


3.61921 


7.99765e-05 


3 


TIT 


6.65201 


0.000489449 


6 


3 x (TT) 


21.4633 


0.000900802 


12 


6 x (TT) 


68.7905 


0.000904738 


20 


io x (TT) 


169.413 


0.00284978 



TABLE I: Selected simulation results, energy and standard deviation 
for different electron configurations and various N- values. 



B. Classical calculations 

A classical mechanics-based "molecular dynamics" ap- 
proach is well suited for massive particles at sufficient high 
temperatures, for example in simulations of molecular liquids. 
Several classical simulations of low temperature and confined 
dilute matter systems have also been successfully applied re- 
cently, see for example Ref . 12111 . We apply the same program 
in the present work except that the ^-direction is kept frozen. 

In essence, the dynamics is described by Newton's equation 
of motion through 



(9) 



where the Force F is derived from Eq. {IJ and modified at 
each time step At according to, 
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The thermostat £ slowly drives the system from a non- 
equilibrium initial configuration to the final state character- 
ized by a required temperature T\ 



(11) 



Here Ek is the kinetic energy, k is Boltzmann's constant and 
7 is a parameter of order unity which plays no role for the fi- 
nal structures, but regulates the computational speed of cool- 
ing/heating of the system. The temperature T' is the set equi- 
librium temperature of the system. 

When propagating the system in time we solve Eq. (|9|) 
numerically with an appropriate time step by the leap-frog 
method 12211 . Time steps and cooling parameters are carefully 
adjusted to avoid the system to be locked into any local energy 
minima close to the absolute minimum as described in 12311 . 



III. RESULTS AND DISCUSSION 
A. Energy levels 

In Fig. [l] we show the resulting ground state energies as 
functions of the number of electrons N and with zero mag- 
netic field and confinement energy Hcoq = 3 meV. The two up- 
per curves are the variational Monte Carlo energies obtained 



without and with the Jastrow factor described in the previ- 
ous section. This factor is seen to play an increasing role for 
increasing number of electrons, as expected since the average 
distance between the electrons decreases with increasing elec- 
trons in the dot. For N = 20 electrons the Jastrow factors de- 
creases the energy with about 10%. To compare with the total 
role of the electron-electron interaction we also plot in Fig. [I] 
the corresponding ground state for non-interacting particles, 
full line with circles. The role of the correlations is again seen 
to increase with the number of electrons, but much less than 
for the interacting case. Finally, in the same figure we display 
the ground state energy obtained with the classical molecular 
dynamics calculations. As expected the energies are a factor 
of two below the quantal interacting energies due to vanishing 
kinetic energy for T = 0. 
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FIG. 1: Ground state energy as a function of number of electrons 
(a = 0, B — 0). Crosses: QMC calculation without Jastrow fac- 
tor; dashed-dotted line: QMC with the Jastrow factor. Full line with 
circles: Independent particle energies. Dashed line: Classical calcu- 
lation. 

All calculations in Fig. [I] are performed with a = 0, see 
Eq. |lj. In Fig.|2|the role of screening is examined by comput- 
ing the addition energy spectra with Hcoq = 3 meV. The addi- 
tion spectra is given by, Sfi = E(N+l)-2E(N)+E(N-i), 
where ji is the chemical potential. The energies have been 
evaluated with a = 0, 5 and 10 nm. Strong peaks at N = 2 
and N = 6 and higher "magic" numbers are obtained in all 
cases, in qualitative agreement with experiment. The strong 
peak for N = 4 gets reduced for increasing a and the calcu- 
lations for a = lOnm are in better agreement with experiment 
than the a = case. One can argue for a finite a value since 
a quantum dot has a finite thickness, resulting in a screened 
Coulomb interaction between the electrons. This leads to a 
reduction in the total energy and makes the harmonic oscilla- 
tor potential more dominant. We also notice that the addition 
energy is large at the half-filled shells N = 4, 9, 12... Accord- 
ing to Hund's rule, these are the states with maximum total 
spin, i.e. the electron eigenspins are parallel. Such configu- 
rations require an antisymmetric spatial wave function, which 
minimises the Coulomb repulsion between the electrons by 



4 



12 



10- \ 




\) 5 10 15 20 

N 



netic fields is of great importance for coherent transport in 
coupled double dots. The tunneling amplitude between dots 
depends exponentially on the barrier parameters and energy 
levels which again can be adiabatically tuned by a time de- 
pendent external field. In Fig. |3 the total energy of a two- 
and one-particle ground state as a function of magnetic field 
strength is shown, in the two-electron case with and without 
the electron-electron interaction. As expected, with zero total 
orbital momentum, the energy increases monotonous with the 
magnetic field. The contribution from the Coulomb repulsion 
is seen to be very significant already for two electrons as it 
accounts for around 40% of the total energy. 

400 | 1 1 1 1 1 1 1 1 1 1 

350 - / - 



FIG. 2: (Color online) Addition energy spectrum, computed for the 
hujQ = 3meV-case, with different offsets a. The blue line corre- 
sponds to a = nm. Points on the green line are computed with 
a = 5 nm, and those marked with x's are with a = 10 nm. The 
points on the black line are those obtained from the experiment of 
Ref. G3|. 



> 

CD 

£ 250 - 
CD 

CD 200 " 
CD 

"CO 150 - 

o 



100 - 



keeping electrons apart and making the configuration more 
stable. By screeening the electrons we reduce the effect of the 
Coulomb repulsion, therefore also making the spatial antisym- 
metrisation less important for lowering the energy of the con- 
figuration. Therefore, the addition energy peaks at half-filled 
shells become smaller as we increase the screening parameter 
a. 
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FIG. 4: Total energy per particle as a function of magnetic field 
strength for N — 6, 12 and 20. The screening parameter is a — 10. 

The importance of the Coulomb repulsion increases with 
increasing magnetic field and number of particles. Introduc- 
ing a finite offset a puts a limit on how close the electrons may 
approach each other. This will of course prevent the Coulomb 
interaction from growing with magnetic field when the offset 
a is large enough. In Fig. [4] we plot the scaled total energy as 
a function of magnetic field strength for N = 6, N = 12 and 
TV = 20 electrons. The comparison shows that the electron 
repulsion amplifies the increase of energies due to the increas- 
ing magnetic field strength. The increased energy is seen to be 
close to proportional with the number of particles. 



B. Electron probability distribution 



FIG. 3: Total energy of a two-particle ground state as a function of 
the strength of the magnetic field. Upper curve: Two particles with 
the electron-electron interaction. Middle curve: Two particles with 
no electron-electron interaction. Lower curve: A single electron. The 
screening parameter is a — 10. 

We now turn to the non-zero magnetic field configurations. 
The detailed lifting of the ground state energies with mag- 



We now turn to a detailed comparison of the quantal ground 
state probability distribution in absence of a magnetic field 
with corresponding classical calculations. In Fig.|5]the classi- 
cal electron configurations (left) with the single-electron prob- 
ability density obtained after integration over the remaining 
N — 1 electron's, p(rjsr) = j dr i... dr n -ll^ n\ 2 (right) are 
shown. The classical frozen state is taken directly from the 
final positions and placed on top of a figure with each elec- 
tron position rotating a number of uniformly distributed an- 
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gles with constant radius. 

For the classical results of Fig.|3 the black dots represent 
electron positions obtained from freezing an initial distribu- 
tion to T' = 0. By keeping T' > 0, but small, we also produce 
a classical spatial density distribution by time averaged super- 
positions of a large number of instantaneous configurations 
as described in Ref . 12^1 . These densities are represented by 
shades of gray superimposed over the static electron positions. 
The TV = 6 configuration is identical to the one obtained in 
Ref. ll23ll . It contains a single electron in the middle and five 
electrons in an outer ring. This double structure is also seen 
for TV = 12, but now with 3 electrons in the inner ring. The 
TV = 19 configuration not shown here but published in l23ll 
contains the first triple ring with one electron in the center, 
six electrons in the middle ring and 12 in the outer ring. It 
is interesting to note that when adding one extra electron, the 
TV = 20 configuration minimizes energy by placing the extra 
electron in the middle ring. 

Comparison between classical and quantal probability den- 
sity distributions shows as expected a clear discrepancy for 
TV = 2 since anti symmetrization results in a spatial symmet- 
ric wave function with a large probability of finding an elec- 
tron in the trap center. For TV = 6 the quantal ground state 
is a single ring with no inner electron, again due to antisym- 
metry, while the classical TV = 6 configuration has a centrally 
positioned electron. This pattern remains valid also for the re- 
maining two cases, viz., where the classical calculations yield 
a centrally placed electron, the quantal calculations show a 
central density dip at the center, while the absence of central 
electron in the classical configuration (TV = 2 and TV = 12) 
is matched by increased central density in the quantal case. In 
the quantal displays the gray shade coding has been selected 
to visualize the density variations which are of the order of 
about 10 - 30% in the TV = 12 and TV = 20 cases. Apart 
from this effect of the Pauli principle the systematics of the 
main shells for TV = 6, TV = 12 and TV = 20 compare well 
in spatial distribution with the corresponding classical ones, 
showing that the classical static configurations may roughly 
describe the ranges of the charge density distributions even 
for systems with relatively few confined electrons. 

The grayscale in quantum densities is such that it reflects 
the density variations which are increasingly smaller and 
smaller with the number of electrons in the inner region (about 
40 percent in the N=6 case, about 20 for n=20 and around 10 
percent for the largest N- values. These variations reflect the 
nodal structure of the last filled oscillator shell. The top or 
dip in the central region follows the parity of the last filled 
shell which is negative for N=6,20 and positive for N= 12,30. 
For N=30 there are 3 visible quantal and classical shells, the 
inner "inversion" is again obvious. Half of the classical elec- 
trons occupy the outer shell somewhat smaller than the quan- 
tal probability density for the reason discussed above. We 
close by concluding that even if there are distinct quantum 
features exposed in the right panels, the shell structure and 
the sizes of the classical systems show a surprisingly good 
correspondence with the quantal structures. 
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FIG. 5: Classical electron configuration (left) and quantum me- 
chanical electron density (right) for closed shell ground states for 
N=2,6,12 20 and 30 electrons. The scale of all figures are from -300 
nm to 300 nm in both directions. 
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IV. CONCLUSION 

The present investigation has demonstrated the usefulness 
of the Quantum Monte Carlo method for studies of relatively 
large numbers electrons inside a quantum dot with parame- 
ters typically achieved in experiments. Studies of the effects 
of the magnetic field have shown that the energy per particle 
increase due to the magnetic field is amplified by the electron 
repulsion, roughly proportional to the number of electrons, his 
follows from comparisons with the non-interacting cases and 
can be traced to the fact that the repulsion leads to increased 
expectation values of the charge distribution radius and the 
fact that the energetic effects of the magnetic field depends on 
this quantity. 

The quantal calculations were compared with classical 
Molecular Dynamics calculations of the static equilibrium po- 



sitions of classical charges confined by a harmonic force. The 
comparisons for closed shells shows that for the chosen type 
of quantum dots, the effect of the Pauli principle leads to dif- 
fering behavior of the classical configurations and the quantal 
distributions in the central area. The overall shell structure 
in the quantal and classical case nevertheless has a remark- 
able agreement in shape and spatial extension. This points 
further to the possibility that a range of important dynami- 
cal processes in quantum dots may be understood from purely 
classical models. 
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